%maximum between-class variance method
%a: image data
function th=thresh_md(a)
%a=imread(p);
%count=a; 
%histogram
a=double(a);
[m,n]=size(a); 
L=max(max(a));
count(1:L+1)=0;
for i=1:n
    count(a(i)+1)=count(a(i)+1)+1;
end
%figure,plot(1:L+1,count);
count=count/n; 
for i=1:L/2+1 
   if count(i)~=0 
     st=i-1; 
     break; 
   end 
end 
for i=L+1:-1:L/2 
   if count(i)~=0 
      nd=i-1; 
      break; 
    end 
end 
%f is the probability of each gray level 
f=count(st+1:nd+1); 
p=st; 
q=nd-st; 
u=0; 
for i=1:q 
   %u is the total mean level 
   u=u+f(i)*(p+i-1); 
   %ua(i) is the first-order cumulative moments of the histogram up to ith
   %level
   ua(i)=u;  
end
u2=0;
for i=1:q
    %u2 is the mean variance level
    u2=u2+(f(i)*(p+i-1)-u).^2;
    u3=0;
    for t=1:i
      u3=u3+(f(t)*(p+t-1)-ua(i)).^2;
    end
    ua2(i)=sqrt(u3/i);
end
%ua
%ua2
u2=sqrt(u2/q);
for i=1:q 
   %w(i)is the cumulative probality of ith pixel
   w(i)=sum(f(1:i));  
end 
%d is the maximum between class variance
d=(u2*w-ua2).^2./(w.*(1-w));
[y,tp]=max(d); 
th=tp+p;
clear p;
clear tp;
clear d;



